In [1]:
# Setting up a model and a mesh for the MT forward problem
In [62]:
import SimPEG as simpeg, sys
import numpy as np
from SimPEG import NSEM
sys.path.append('/Volumes/MacintoshHD/Users/thibautastic/Dropbox/PhD_UBC/telluricpy')
import telluricpy
import matplotlib.pyplot as plt
%matplotlib inline
import copy
In [3]:
# Define the area of interest
bw, be = 557100, 557580
bs, bn = 7133340, 7133960
bb, bt = 0,480
hSize,vSize = 25., 10 nrCcore = [15, 8, 6, 5, 4, 2, 2, 2, 2] hPad = simpeg.Utils.meshTensor([(hSize,10,1.5)]) hx = np.concatenate((hPad[::-1],np.ones(((be-bw)/hSize,))hSize,hPad)) hy = np.concatenate((hPad[::-1],np.ones(((bn-bs)/hSize,))hSize,hPad)) airPad = simpeg.Utils.meshTensor([(vSize,13,1.5)]) vCore = np.concatenate([ np.ones(i)*s for i, s in zip(nrCcore,(simpeg.Utils.meshTensor([(vSize,1),(vSize,8,1.3)])))])[::-1] botPad = simpeg.Utils.meshTensor([(vCore[0],8,-1.5)]) hz = np.concatenate((botPad,vCore,airPad))
x0 = np.array([bw-np.sum(hPad),bs-np.sum(hPad),bt-np.sum(vCore)-np.sum(botPad)])
meshFor = simpeg.Mesh.TensorMesh([hx,hy,hz],x0)
In [76]:
# Build the Inversion mesh
# Design the tensors
hSizeI,vSizeI = 25., 10.
nrCcoreI = [12, 6, 6, 6, 5, 4, 3, 2, 1]
hPadI = simpeg.Utils.meshTensor([(hSizeI,9,1.5)])
hxI = np.concatenate((hPadI[::-1],np.ones(((be-bw)/hSizeI,))*hSizeI,hPadI))
hyI = np.concatenate((hPadI[::-1],np.ones(((bn-bs)/hSizeI,))*hSizeI,hPadI))
airPadI = simpeg.Utils.meshTensor([(vSizeI,12,1.5)])
vCoreI = np.concatenate([ np.ones(i)*s for i, s in zip(nrCcoreI,(simpeg.Utils.meshTensor([(vSizeI,1),(vSizeI,8,1.3)])))])[::-1]
botPadI = simpeg.Utils.meshTensor([(vCoreI[0],7,-1.5)])
hzI = np.concatenate((botPadI,vCoreI,airPadI))
# Calculate the x0 point
x0I = np.array([bw-np.sum(hPadI),bs-np.sum(hPadI),bt-np.sum(vCoreI)-np.sum(botPadI)])
# Make the mesh
meshInv = simpeg.Mesh.TensorMesh([hxI,hyI,hzI],x0I)
meshFor = copy.deepcopy(meshInv)
In [77]:
print np.sum(vCore)
print meshFor.nC
print meshFor
In [78]:
# Save the mesh
meshFor.writeVTK('nsmesh_CoarseHKPK1.vtr',{'id':np.arange(meshFor.nC)})
nsvtr = telluricpy.vtkTools.io.readVTRFile('nsmesh_CoarseHKPK1.vtr')
In [79]:
nsvtr
Out[79]:
In [80]:
topoSurf = telluricpy.vtkTools.polydata.normFilter(telluricpy.vtkTools.io.readVTPFile('../../Geological_model/CDED_Lake_Coarse.vtp'))
activeMod = telluricpy.vtkTools.extraction.extractDataSetWithPolygon(nsvtr,topoSurf)
In [81]:
topoSurf
Out[81]:
In [82]:
#telluricpy.vtkTools.io.writeVTUFile('activeModel.vtu',activeMod)
In [83]:
# Get active indieces
activeInd = telluricpy.vtkTools.dataset.getDataArray(activeMod,'id')
In [84]:
# Make the conductivity dictionary
# Note: using the background value for the till, since the extraction gets the ind's below the till surface
geoStructFileDict = {'Till':1e-4,
'PK1':5e-2,
'HK1':1e-3,
'VK':5e-3}
In [85]:
# Loop through
extP = '../../Geological_model/'
geoStructIndDict = {}
for key, val in geoStructFileDict.iteritems():
geoPoly = telluricpy.vtkTools.polydata.normFilter(telluricpy.vtkTools.io.readVTPFile(extP+key+'.vtp'))
modStruct = telluricpy.vtkTools.extraction.extractDataSetWithPolygon(activeMod,geoPoly,extBoundaryCells=True,extInside=True,extractBounds=True)
geoStructIndDict[key] = telluricpy.vtkTools.dataset.getDataArray(modStruct,'id')
In [86]:
# Make the physical prop
sigma = np.ones(meshFor.nC)*1e-8
sigma[activeInd] = 1e-3 # 1e-4 is the background and 1e-3 is the till value
# Add the structure
for key in ['Till','PK1','HK1','VK']:
sigma[geoStructIndDict[key]] = geoStructFileDict[key]
In [87]:
fig = plt.figure(figsize=(10,10))
ax = plt.subplot(111)
model = sigma.reshape(meshFor.vnC,order='F')
a = ax.pcolormesh(meshFor.gridCC[:,0].reshape(meshFor.vnC,order='F')[:,20,:],meshFor.gridCC[:,2].reshape(meshFor.vnC,order='F')[:,20,:],np.log10(model[:,20,:]),edgecolor='k')
ax.set_xlim([bw,be])
ax.set_ylim([-0,bt])
ax.grid(which="major")
plt.colorbar(a)
ax.set_aspect("equal")
In [88]:
# Save the model
meshFor.writeVTK('nsmesh_CoarseHKPK1_NoExtension.vtr',{'S/m':sigma})
In [89]:
import numpy as np
In [90]:
# Set up the forward modeling
freq = np.logspace(5,3,16)
np.save('MTfrequencies',freq)
In [91]:
freq
Out[91]:
In [ ]:
In [92]:
# Find the locations on the surface of the model.
# Get the outer shell of the model
import vtk
actModVTP = telluricpy.vtkTools.polydata.normFilter(telluricpy.vtkTools.extraction.geometryFilt(activeMod))
polyBox = vtk.vtkCubeSource()
polyBox.SetBounds(bw-5.,be+5,bs-5.,bn+5.,bb-5.,bt+5)
polyBox.Update()
# Exract the topo of the model
modTopoVTP = telluricpy.vtkTools.extraction.extractDataSetWithPolygon(actModVTP,telluricpy.vtkTools.polydata.normFilter(polyBox.GetOutput()),extractBounds=False)
telluricpy.vtkTools.io.writeVTPFile('topoSurf.vtp',actModVTP)
In [93]:
# Make the rxLocations file
x,y = np.meshgrid(np.arange(bw+12.5,be,25),np.arange(bs+12.5,bn,25))
xy = np.hstack((x.reshape(-1,1),y.reshape(-1,1)))
# Find the location array
locArr = telluricpy.modelTools.surfaceIntersect.findZofXYOnPolydata(xy,actModVTP) #modTopoVTP)
np.save('MTlocations',locArr)
In [94]:
telluricpy.vtkTools.io.writeVTPFile('MTloc.vtp',telluricpy.dataFiles.XYZtools.makeCylinderPtsVTP(locArr,5,10,10))
In [95]:
# Running the forward modelling on the Cluster.
# Define the forward run in findDiam_MTforward.py
In [57]:
%matplotlib qt
#sys.path.append('/home/gudni/Dropbox/code/python/MTview/')
import interactivePlotFunctions as iPf
In [60]:
# Load the data
mtData = np.load('MTdataStArr_nsmesh_HKPK1Coarse_noExtension.npy')
mtData
Out[60]:
In [61]:
In [110]:
iPf.MTinteractiveMap([mtData])
Out[110]:
In [104]:
# Looking at the data shows that data below 100Hz is affected by the boundary conditions,
# which makes sense for very conductive conditions as we have.
# Invert data in the 1e5-1e2 range.
Run the inversion on the cluster using the inv3d/run1/findDiam_inversion.py
In [106]:
drecAll = np.load('MTdataStArr_nsmesh_0.npy')
In [109]:
np.unique(drecAll['freq'])[10::]
Out[109]:
In [93]:
# Build the Inversion mesh
# Design the tensors
hSizeI,vSizeI = 25., 10.
nrCcoreI = [12, 6, 6, 6, 5, 4, 3, 2, 1]
hPadI = simpeg.Utils.meshTensor([(hSizeI,9,1.5)])
hxI = np.concatenate((hPadI[::-1],np.ones(((be-bw)/hSizeI,))*hSizeI,hPadI))
hyI = np.concatenate((hPadI[::-1],np.ones(((bn-bs)/hSizeI,))*hSizeI,hPadI))
airPadI = simpeg.Utils.meshTensor([(vSizeI,12,1.5)])
vCoreI = np.concatenate([ np.ones(i)*s for i, s in zip(nrCcoreI,(simpeg.Utils.meshTensor([(vSizeI,1),(vSizeI,8,1.3)])))])[::-1]
botPadI = simpeg.Utils.meshTensor([(vCoreI[0],7,-1.5)])
hzI = np.concatenate((botPadI,vCoreI,airPadI))
# Calculate the x0 point
x0I = np.array([bw-np.sum(hPadI),bs-np.sum(hPadI),bt-np.sum(vCoreI)-np.sum(botPadI)])
# Make the mesh
meshInv = simpeg.Mesh.TensorMesh([hxI,hyI,hzI],x0I)
In [94]:
meshInv.writeVTK('nsmesh_HPVK1_inv.vtr',{'id':np.arange(meshInv.nC)})
In [101]:
nsInvvtr = telluricpy.vtkTools.io.readVTRFile('nsmesh_HPVK1_inv.vtr')
activeModInv = telluricpy.vtkTools.extraction.extractDataSetWithPolygon(nsInvvtr,topoSurf,extBoundaryCells=True)
In [102]:
sigma = np.ones(meshInv.nC)*1e-8
indAct = telluricpy.vtkTools.dataset.getDataArray(activeModInv,'id')
sigma[indAct] = 1e-4
In [103]:
meshInv.writeVTK('nsmesh_HPVK1_inv.vtr',{'id':np.arange(meshInv.nC),'S/m':sigma})
In [108]:
from pymatsolver import MumpsSolver
In [106]:
pymatsolver.AvailableSolvers
Out[106]:
In [116]:
NSEM.Utils.skindepth(1000,100000)
Out[116]:
In [117]:
np.unique(mtData['freq'])[5::]
Out[117]:
In [ ]: